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ABSTRACT 

We have performed a large set of high-resolution cosmological simulations using 
smoothed particle hydrodynamics (SPH) to study the formation of the first luminous 
objects in the ACDM cosmology. We follow the collapse of primordial gas clouds in 
eight early structures and document the scatter in the properties of the first star- 
forming clouds. Our first objects span formation redshifts from z ~ 10 to z ~ 50 and 
cover an order of magnitude in halo mass. We find that the physical properties of 
the central star- forming clouds are very similar in all of the simulated objects despite 
significant differences in formation redshift and environment. This suggests that the 
formation path of the first stars is largely independent of the collapse redshift; the 
physical properties of the clouds have little correlation with spin, mass, or assembly 
history of the host halo. The collapse of proto-stellar objects at higher redshifts pro- 
gresses much more rapidly due to the higher densities, which accelerates the formation 
of molecular hydrogen, enhances initial cooling and shortens the dynamical timescales. 
The mass of the star- forming clouds cover a broad range, from a few hundred to a 
few thousand solar masses, and exhibit various morphologies: some have disk-like 
structures which are nearly rotational supported; others form flattened spheroids; still 
others form bars. All of them develop a single protostellar 'seed' which does not frag- 
ment into multiple objects up to the moment that the central gas becomes optically 
thick to H2 cooling lines. At this time, the instantaneous mass accretion rate onto the 
centre varies significantly from object to object, with disk-like structures having the 
smallest mass accretion rates. The formation epoch and properties of the star-forming 
clouds are sensitive to the values of cosmological parameters. 

Key words: methods: cosmology: theory - early universe - galaxies: formation - 
hydrodynamics - molecule processes - stars: formation 



1 INTRODUCTION 

The first generation stars, often referred to as "Population- 
Ill" stars, are thought to be the first sources of light in our 
Universe and the origin of the heavy elements required for 
the subsequent formation of ordinary stellar populations. 
Their remnant black holes may be the seeds from which 
supermassive black holes grow, including those that power 
quasars at redshift ~ 6 (e.g. Fan et al. 2003), when the 
Universe was just ~ 800 million years old. In the standard 
ACDM cosmology with a high amplitude of initial density 
perturbations, as = 0.9, the first stars are predicted to form 
at redshifts of about 20 and higher, based on the argument 
that gas in 3cr-halos at such redshifts should be able to cool 
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and condense by molecular hydrogen cooling. Recent deter- 
minations of the cosmological parameters, obtained by com- 
bining the 2dFGRS galaxy power spectrum with the first 
and/or third year WMAP data (Sanchez et al 2005, Spergel 
et al. 2006), give a lower value of as ~ 0.75, implying that 
the formation epoch of the first stars could be as late as 
z ~ 10, opening the possibility that these stars may be de- 
tected directly by the next generation of space-borne and 
ground-based telescopes. 

From a theoretical point of view, the formation of the 
first stars is a well-defined problem. Quantum fluctuations 
imprinted during the inflationary epoch provide the initial 
seeds for the formation and growth of the dark matter halos 
that host the first stars. In the absence of heavy elements, 
the only efficient coolant that can promote the formation 
of primordial star-forming gas clouds is molecular hydro- 
gen (Peebles & Dicke 1968; Matsuda, Sato & Takeda 1969). 
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Thus once the cosmological parameters and the fluctuation 
spectrum are specified, the formation of the first luminous 
objects in the universe becomes, in principle, a straightfor- 
ward physics problem. 

Over the past decades, much progress has been achieved 
in following this program. Early studies used semi-analytical 
methods (e.g. Doroshkevich, Zel'Dovich & Novikov 1967; 
Peebles & Dicke 1968; Couchman & Rees 1986; Tegmark 
et al. 1997; Silk & Langer 2006), while others utilised di- 
rect numerical simulations, either in one dimension (Haiman 
& Thoul & Loeb 1996; Omukai & Nishi 1998; Nakamura 
& Umemura 2001; Ahn & Shapiro 2007) or in three di- 
mensions (Abel, Bryan & Norman 2000, 2002, hereafter 
ABN02; Bromm et al. 1999, 2002; Fuller & Couchman 2000; 
Yoshida et al. 2003, 2006; O'Shea & Norman 2006a,b). While 
the semi-analytical calculations have been instrumental in 
sketching out the relevant physics, the numerical simula- 
tions have been able to deliver detailed predictions for the 
outcome of the coupling between the chemo-thermal evolu- 
tion and gravitational collapse. Comprehensive reviews of 
recent progress on the formation of the first luminous ob- 
jects and associated radiative feedback processes may be 
found in Barkana & Loeb (2001), Bromm & Larson (2004), 
Glover (2005) and Ciardi & Ferrara (2005), and Ripamonti 
& Abel (2005). 

Pioneering three-dimensional simulations of the forma- 
tion of the first objects including non-equilibrium chemistry 
suggested that the first stars could be more massive than 
100 M (Abel, Bran & Norman 2000; 2001; Bromm et al. 
1999, 2002;). ABN02 employed cosmological initial condi- 
tions and adaptive mesh refinement (AMR) techniques to 
follow a star-forming gas cloud until it reached a central hy- 
drogen number density, nu ~ 10 10 cm -3 , the point at which 
the assumption that the gas is everywhere optically thin 
breaks down and a radiative transfer calculation is required. 
Recently, Yoshida et al. (2006) included the effect of molec- 
ular line opacities and followed the further collapse of the 
central core to nn ~ 10 16 cm -3 . Bromm et al. (2002) (here- 
after BCL02) and Bromm & Loeb (2004) used constrained 
initial conditions and smooth particle hydrodynamics (SPH) 
to study the same problem. Their simulations agree with 
those of ABN02 on many of the properties of the result- 
ing star- forming cloud, but they also differ in a number of 
details. 

A major disagreement is that ABN02 proposed that 
only one star forms per halo, whilst BCL02 argued that 
multiple stellar systems can form if the gas has sufficient 
initial rotation. In the latter case, the gas first collapses into 
a disk-like structure which fragments into several clumps 
very quickly. These clumps will then make their own stars 
if they can survive negative feedback effects from more ad- 
vanced neighbouring clumps. Since ABN02 only simulated 
one particular object and the initial set-up of BCL was not 
realistic, there is hope that progress can be made by carrying 
out a number of simulations with realistic initial conditions 
in order to establish the typical properties of star-forming 
clouds and their variation. 

One important caveat for simulations of the first stars 
has been noted by White & Springel (2000) and Gao et 
al. (2005) (hereafter G05). They argue that the first ha- 
los of any given mass do not form in 'typical' regions, but 
rather in 'protocluster' regions corresponding to peaks of 



the large-scale density field in the early universe. Since a 
significant contribution to the local overdensity at the sites 
of first star formation comes from megaparsec and larger 
wavelength fluctuations, some conclusions about the first 
objects reached from simulations of small volumes (which 
represent 'typical' regions) may be misleading. In particu- 
lar, the formation redshift of the first halos of any given mass 
or temperature will be systematically underestimated. 

G05 carried out a sequence of Af-body re- simulations 
to follow the growth of one of the most massive progeni- 
tors of a supercluster region, from redshift z ~ 80, when 
its mass was just 10 M©, to the present epoch. They found 
that the mass of the first object reaches 1O 5 M0 already at 
z ~ 50. A follow-up study of Reed et al. (2005) grafted a 
semi-analytical model developed by Yoshida et al. (2003) 
onto the dark matter merger trees of the G05 simulations 
and concluded that the first halo in this sequence of sim- 
ulations would make one or more stars already at z ~ 47. 
Very rare objects like the one simulated by G05 are differ- 
ent from those formed in typical regions. These rare objects 
reside in extremely dense environments and experience par- 
ticularly rapid growth. Their higher collapse redshift results 
in a higher initial density than in a typical region and this 
could result in significantly different evolution later on be- 
cause the dynamical time, the cooling time, and the reaction 
rates for the formation of molecular hydrogen all depend on 
gas density. It is therefore an interesting question whether 
the formation paths and properties of star-forming clouds 
in these very rare objects are at all similar to those forming 
later on in more typical environments. 

In this paper, we address this issue by asking two impor- 
tant questions: (1) Is there a significant scatter in the prop- 
erties of primordial star-forming clouds? (2) Do the prop- 
erties of star-forming clouds depend on formation redshift? 
We perform simulations of a sample of 8 first objects. By fol- 
lowing evolution from different sets of initial conditions we 
can simulate a variety of formation redshifts and investigate 
the redshift dependence of the formation process. 

This paper is structured as follows. In Section 2, we 
briefly introduce the chemical processes that are important 
for primordial gas. In Section 3, we describe our simulations. 
In Section 4, we take three of them as examples to sketch a 
general picture for the formation path of the first stars and 
then focus on one of them in greater detail. In Section 5, 
we study whether the properties of the star-forming clouds 
depend on collapse redshift and, in Section 6, we investigate 
how the baryon fraction influences the formation of the first 
stars. In Section 7, we discuss the abundance of star-forming 
halos. Finally, we give a summary of our results and set out 
our conclusions in Section 8. 



2 PRISTINE GAS CHEMISTRY 

In this section, we briefly summarize the basic gas processes 
that are important for primordial gas unpolluted by metals. 
Full details are given by Abel et al. (1997), Galli & Palla 
(1998) and Yoshida et al. (2006). For simulations that start 
at very high z, we now also include photo-reactions with 
cosmic microwave background (CMB) photons. At high red- 
shift, in the absence of metals and when the temperature of 
gas in halos is lower than 10 4 K (above which atomic hydro- 
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gen line cooling is dominant), the principal coolant that can 
enable gas to condense into stars is the H 2 molecule. 

The formation of molecular hydrogen has 3 main chan- 
nels. At high redshifts, z > 200, the channel given by 



: + + h 

i + h 



H 2 + + 7 
H 2 + H+ 



is dominant (but see Hirata & Padmanabhan (2006) for a 
different viewpoint). At these early times, there are suffi- 
cient numbers of energetic CMB photons to destroy H~ ions 
rapidly, so that a second possible channel, 



H + e" 
H" + H 



H" + 7 
H 2 + e~ 



is initially strongly suppressed. However, as the Universe ex- 
pands and the CMB photons become less energetic, the H~ 
channel eventually becomes the main mechanism for creat- 
ing molecular hydrogen. In this process, the residual elec- 
trons left unattached after recombination act as a catalyst. 

When the density of gas exceeds nu ~ 10 8 cm -3 , a third 
channel becomes important, through the three-body reac- 
tion (Palla, Salpeter & Stahler 1983): 

H + H + H ^ H 2 + H, 
H + H + H 2 -> H 2 + H 2 . 

At sufficiently high densities, this reaction can, in fact, con- 
vert almost all atomic hydrogen into molecular form in a 
very short time. Our simulation code follows reactions for 9 
species (e~, H, H + , He, He + , He ++ , H 2 , Hj, H~), including 
these H2 formation processes. 

The cooling function for molecular hydrogen can be ex- 
pressed as 



A(n H ,T) 



A L te(T) 



1 + 



Alte(T) 



(1) 



n H A n 



->0(T) 



where Alte and A nH ^o are the cooling functions for high 
and low density limits. We adopt the high density limit 
given by Hollenbach & McKee (1979), assuming local ther- 
mal equilibrium (LTE), and the low density limit of Galli 
& Palla (1998). The transition between these two regimes 
occurs approximately at a gas density of nu ~ 10 4 cm -3 . 
Above this value, the cooling rate per hydrogen molecule is 
independent of density, while below it, it is proportional to 
density. 



3 SIMULATIONS 

We have carried out a suite of simulations designed to study 
the formation of the first stars. For the first series, a stan- 
dard ACDM cosmology, with matter, baryon and dark en- 
ergy density parameters, Q m = 0.3, = 0.04, Qa = 0.7 re- 
spectively, and a Hubble constant /iioo = 0.7, was adopted. 
The initial power spectrum was generated with CMBFAST 
and normalised to give an rms linear fluctuation amplitude 
of as = 0.9 in 8 /i _1 Mpc tophat spheres at z = 0. The initial 
ionisation fraction was computed using RECFAST (Seager, 
Sasselov & Scott 2000). 

One of the simulations used in this study has been pre- 
viously analysed extensively by G05 and Reed et al. (2005). 



This simulation was constructed through a sequence of re- 
simulations of the most massive progenitors of a superclus- 
ter region. Here we briefly describe the multi-grid procedure 
devised by G05: 

(1) Identify a rich cluster size halo at z = in a cosmolog- 
ical simulation of a very large volume (~ lGpc 3 ). 

(2) Resimulate the evolution of this rich cluster and its 
environment with higher mass resolution. 

(3) Determine the most massive object in the high reso- 
lution region of the re-simulation as a function of redshift, 
and identify the time when it first contains more than 10000 
particles inside the viral radius. 

(4) Resimulate the evolution of this progenitor object and 
its immediate surroundings with a further improvement in 
mass and force resolution. 

(5) Iterate steps 3 and 4 until the desired redshift and pro- 
genitor mass are reached. 

Following this procedure through 5 levels of re- 
simulation, G05 were able to resolve a massive progenitor 
of the supercluster region capable of efficient H 2 -molecular 
cooling at redshift z ~ 50. Their final, highest-resolution 
simulation had a particle mass of 1.08M© and is labelled 
"R5." Here we have re-ran "R5," this time including ra- 
diative cooling by molecular hydrogen and non-equilibrium 
chemistry. This is one of the main simulations analysed in 
this paper. 

We carried out two additional simulations drawn from a 
periodic box of 300 /i _1 kpc on a side. As noted by G05, the 
use of periodic boundary conditions suppresses fluctuations 
with wavelength longer than the side of the computational 
box. Therefore, if the box size is small, the appearance of a 
halo that could host first stars will be significantly delayed 
compared to a simulation of a larger volume. Our main mo- 
tivation for simulating these two small periodic boxes is to 
enable us to compare our results to earlier studies carried 
out in this way (e.g. ABN02). In addition, however, since the 
formation of the first stars could plausibly span a wide range 
in redshift, our two small periodic box simulations are still 
interesting for studying first star formation at comparatively 
low redshifts. 

In order to achieve sufficient resolution, we again 
utilised a resimulation technique similar to that of G05. 
First, we performed a simulation of a 300/i _1 kpc cubic re- 
gion using 256 3 particles. We then identified the most mas- 
sive objects at redshift z ~ 15 which typically have mass 
of ~ 1O 6 M and resimulated them with a resolution com- 
parable to the R5 simulation. Note that since the dynamic 
range of this simulation is much smaller than that of R5 (the 
initial volume is a factor 4 x 10 9 smaller), only one level of 
refinement, rather than five, is required to create the initial 
conditions of the final target simulation. Table 1 summarizes 
the principal numerical parameters of these simulations. 

Since uncertaintes remain on the values of the cosmolog- 
ical parameters which describe our Universe, we have also re- 
run all the simulations listed in Table 1 adopting the param- 
eters given in Table 5 of Spergel et al. (2006), derived from 
the combination of 3 years WMAP microwave data with the 
2dFGRS galaxy power spectrum (Cole et al. 2005). We re- 
tain the same phases in the initial conditions but adjust the 
amplitudes of the modes according to the new power spec- 
trum shape and normalisation. These cosmological param- 
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R5 Zl Z2 



^dm[M ] 


1.08 


3.1 


2.12 


m g as [M©] 


0.166 


0.48 


0.326 


M 2O o[M ] 


2.21 x 10 5 


3.75 x 10 5 


6.18 x 10 5 


^end 


47.91 


25.9 


21.0 



Table 1. Numerical parameters of our first set of simulations. 
Here md m is the mass of each dark matter particle, while m gas is 
the mass of each gas particle; z e nd is the time when the central 
density of the gas cloud reaches ~ 10 10 when the optically thin 
assumption breaks down; M200 is the halo mass defined as the 
enclosed mass in a sphere of mean interior overdensity 200 times 
the cosmic mean density at the final time. This simulation se- 
ries was performed using the standard cosmological parameters, 
Qm = 0.3, Q A = 0.7, Q b = 0.04, h 100 = 0.7, a 8 = 0.9, with a 
primordial power spectrum index of n s = 1. 





R5wt 


Zlwt 


Z2wt 


ra dm [M©] 


0.77 


2.22 


1.52 


ragas[M ] 


0.167 


0.48 


0.328 


M 2O o[M ] 


4.42 x 10 5 


1.16 x 10 6 


2.55 x 10 6 


z end 


30.04 


13.43 


10.33 



Table 2. Numerical parameters of our second set of simula- 
tions. The values of the cosmological parameters were taken from 
Table 5 of Spergel et al. (2006): ft m = 0.234, ^ A = 0.764, 
tt b = 0.042, h 100 = 0.73, cr 8 = 0.74, and n s = 0.95. The meaning 
of the tabulated quantities is the same as in Table 1. 

eters are significantly different from the values adopted in 
our other simulations: Q m = 0.236, Q A = 0.764, Q b = 0.042 
and h±oo = 0.73. The primordial power spectrum slope in- 
dex is n s = 0.95, and the normalisation is as = 0.74. Further 
details of this set of simulations are given in Table 2. 

We performed a third set of controlled simulations to 
complement the other two. By design, these simulations 
form stars at redshifts intermediate between those of the 
first and the second sets. To achieve this, we ran additional 
versions of Zl and Z2, with a lower value of as = 0.74, but 
employing the standard power spectrum, in order to increase 
the power on small scales, whilst the other cosmological pa- 
rameters were set according to the 3rd year WMAP results. 
These simulations are detailed in Table 3. 

Finally, we carried out simulations with both a higher 
and a lower baryon fraction than in the Zl simulation in 
order to investigate how this parameter influences the for- 





Zlw 


Z2w 


m dm [M Q ] 


2.22 


1.52 


m g as [M ] 


0.48 


0.328 


M 2OO [M ] 


4.94 x 10 5 


6.9 x 10 5 


z end 


20.65 


16.75 



Table 3. Parameters of simulations where the three- year WMAP 
cosmological parameters are adopted, but the shape of a standard 
cosmology power spectrum is retained. 





Zlhb 


Zllb 


Wdm[M ] 


2.94 


3.36 


m g as [M©] 


0.63 


0.21 


M 200 [M ] 


3.25 x 10 5 


6.9 x 10 5 


fb 


0.178 


0.06 


z end 


26.26 


24.26 



Table 4. Numerical parameters of our simulations where we ex- 
plored a higher and lower baryon fraction version than used in 
our default simulation of the Zl object. The standard cosmology 
is assumed. 

mat ion of the first stars. For the higher baryon fraction case, 
we chose f b = Q b /Q m — 0.178 which is consistent with the 
3-year WMAP cosmology while, for the lower baryon frac- 
tion case, we adopted f b = 0.06 which matches the value 
used by ABN02 (for their assumed standard Q m = 1 CDM 
cosmology). This allows a more direct comparison of their 
results with ours. See Table 4 for more information about 
these two simulations. 

For all our simulations we used a fixed comoving grav- 
itational softening length for dark matter particles equal to 
1/30 of the mean interparticle separation, while for gas par- 
ticles the gravitational softening was set equal to the SPH 
smoothing length. We experimented with a fixed gravita- 
tional softening length for gas particles, but this caused 
serious numerical problems for collapse simulations such 
as those presented in this paper (Bate & Burkert 1997; 
Abel et al., in preparation). Our simulations were performed 
with the widely used cosmological TreeSPH-code GADGET2 
(Springel 2005). 



4 FORMATION PATH OF THE FIRST STARS 

After recombination, the temperature of the cosmic gas re- 
mains initially coupled to that of the radiation due to the 
presence of a small residual fraction of free electrons. Af- 
ter z rsj 120, this coupling becomes weak and the inter- 
galactic medium (IGM) temperature starts to fall adiabati- 
cally as T oc (1 + z) 2 (e.g. Peebles 1993). While the assem- 
bly of cold dark matter into halos proceeds hierarchically 
starting at tiny mass scales, the accretion of gas into ha- 
los is affected by thermal pressure (e.g. Yoshida, Sugiyama, 
& Hernquist 2003). The density of gas only rises signifi- 
cantly above the mean cosmic gas density within halos that 
have a mass comparable or greater than the Jeans' mass 
Mj = 4/37r 5/2 ciJ(G'p)~ 1/2 , where c s is the sound speed of 
the gas, and p is the density. 

Once a dark matter halo becomes massive enough for 
its baryon content to overcome the thermal gas pressure, it 
can accrete enough material to reach a baryon fraction close 
to the universal value. When the temperature is high enough 
(rsj 1000K), the processes described in Section 2 can promote 
the formation of molecular hydrogen which then acts as a 
coolant. A runaway cooling and compression process begins 
at the halo centre, the details of which are the main topic 
of this and the next section. When the central density of 
hydrogen reaches values of around nu ~ 10 10 cm -3 , the gas 
ceases to be optically thin to the molecular hydrogen lines. 
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Figure 1. Evolution of different physical properties of the 100 SPH particles (~ 20 M©) which end up closest to the centre of the halo 
at the final output time. Panel (a): the black solid lines with symbols show the virial mass of the main halo progenitor; the lower (red) 
solid lines show the mass of the gas component within the virial radius of the halo; the (red) dashed lines are the fiducial gas mass that 
the halo would contain at the universal baryon fraction, fa = Q{,/f2 m . Panel (b): the evolution of the H2-fraction for the 100 selected gas 
particles. Panel (c): the gas temperature evolution. The lower (blue) solid line is an estimate of the virial temperature of the main halo, 
as described in the text; the black solid lines are the maximum temperature defined as the maximum of the spherically mass- averaged 
temperature profile of the the gas in the main halo. The (red) solid lines with symbols show the evolution of the median value of the 
temperature for the 100 selected particles. Panel (d): the solid lines with symbols are the median values of the number density of hydrogen 
nuclei in the 100 selected particles. The dashed lines indicate 200 times the cosmological mean number density. 



At this point, our simulations are no longer able to follow 
the collapse process any further in a physically meaningful 
way because the appropriate radiative transfer processes are 
not modelled. For the purposes of our analysis we therefore 
define the final time of our simulations to be the point where 
the central density reaches a value of nu ~ 10 10 cm -3 . 

In the following subsections, we illustrate the path lead- 
ing to the first stars by discussing three of our simulations, 
R5, Zl, and Z2wt, in detail. These form a comprehensive 
set in terms of the formation redshift of their first star: R5 
is the first to produce a star- forming object, at a redshift of 
z ~ 50; in Zl the halo collapses at an intermediate redshift 
of z ~ 26, close to what is often assumed in the literature to 
be the typical formation epoch of the first stars; Z2wt forms 
its first star as late as z ~ 10. 



4.1 Overview of the formation of the first stars at 
the highest redshifts 

A useful way to describe the formation path of the first stars 
is to track the evolution of the properties of gas particles that 
end up in the central region at the final time. To this end, 
we select the 100 innermost SPH particles which lie closest 
to the centre of the R5 halo at the end of the simulation. We 
take the centre of a halo to be the position of the minimum 
of the gravitational potential. Our results are quite similar 
when we consider the innermost 1000 particles instead. 

The median values of various physical properties of the 



innermost 100 particles as a function of time are shown in 
Figure 1. Panel (a) shows the mass accretion history of the 
gas and dark matter components of the object in R5. The 
solid black line shows the evolution of the virial mass of the 
dominant halo. The virial mass is defined as the total mass 
within the spherical region around the centre that encloses a 
density of 200 times the mean cosmic value. The red dashed 
line is a fiducial gas mass, given by the gas mass that the 
dark matter halo would contain if it had the universal baryon 
fraction, /& = Qb/^m- The solid red line shows the actual 
gas mass within the virial radius. It is not surprising that at 
z > 55 the baryon fraction of the halo is much lower than the 
cosmic mean, because the gas thermal pressure keeps the gas 
out of the halo. The baryon fraction is less than 10% of the 
fiducial value at redshift z — 60 when the halo has a mass of 
1000 M , but rises to 65% at redshift z = 50 by which time 
the mass of the halo has grown to 10000 M©. The universal 
baryon fraction is approximately reached when the halo has 
grown to a mass of 10 5 M©, at z ~ 49. 

Panel (b) shows the evolution of the molecular hydrogen 
fraction, /h 25 of the 100 selected gas particles, while panel 
(c) shows their temperature. For comparison, the blue solid 
line shows the estimated virial temperature of the halo, 

T = nm p V c 2 /(2k B ), (2) 

where \i m p is the mean molecular weight of the gas, V c is the 
circular velocity of the halo and fcs is Boltzmann's constant. 
The black solid line refers to the maximum temperature of 
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the main object, determined by first binning the gas parti- 
cles within the halo into concentric logarithmically spaced 
shells centred on the potential minimum, then computing 
the mass- weighted temperature of each bin, and finally tak- 
ing the maximum value of all the bins. The location of the 
maximum temperature roughly indicates where the accre- 
tion shock occurs. 

The evolution of the number density of hydrogen nu- 
clei for the selected 100 particles is shown in panel (d) as a 
solid line with symbols. The dashed line corresponds to 200 
times the cosmic mean number density. The crossing point 
of the two curves roughly marks the time when these 100 
gas particles settle into the virialised halo. 

Based on these plots, we can outline a simple picture 
of how the gas cloud condenses to the point of forming the 
first star in the R5 halo. Prior to redshift z ~ 53, when the 
mass of the dark halo is less than 2 x 10 4 M©, the virial 
temperature of the halo is substantially lower than that of 
the ambient IGM gas. (Note that the mean temperature 
of the surrounding IGM is higher than the cosmic mean of 
T rsj 60 K because the R5 object resides in a significantly 
overdense region; see G05). In this phase, the gas accretion 
onto the R5 halo is limited by gas pressure. At z ~ 53, a 
major merger takes place and the mass of the halo nearly 
doubles. This significant increase in mass deepens the grav- 
itational potential well, allowing gravity to overwhelm the 
gas pressure, so that ambient gas starts to settle into the 
halo. This can be observed in panel (a) as a rapid increase 
in the baryon fraction at this redshift. The 100 selected gas 
particles settle into the halo as their density reaches the 
characteristic value of 200 times the cosmic mean required 
for virialization, as seen in panel (d). 

The amount of molecular hydrogen increases by almost 
one order of magnitude from z = 53 to z = 50 due to the 
substantial increase in gas temperature as well as density. 
Further production of molecular hydrogen, however, is re- 
tarded by the depletion of electrons due to recombination, 
but the gas nevertheless continues to contract. Shortly be- 
fore the gas density reaches nu ~ 10 8 cm -3 , three-body for- 
mation of molecular hydrogen becomes important and this 
results in very efficient cooling which exacerbates the run- 
away collapse. When the gas density reaches ~ 10 10 cm -3 , 
the gas becomes optically thick to H2-cooling lines, and our 
simulation cannot proceed further in a physically meaningful 
way (see Yoshida et al. 2006 for a calculation of gas evolution 
beyond this point). 

4.2 The formation of first stars at lower redshifts 

(z = 25 and z 10) 

The simulations labelled Zl and Z2wt in Section 3 allow us 
to examine whether or not the formation paths of protostars 
at lower redshifts differ from those forming at high z. The 
first objects in these simulations form much later than in 
R5. In Zl the first protostar forms at z ~ 25, which is often 
taken in the literature to be the typical formation redshift of 
primordial stars; in Zlwt the first protostar forms at z ~ 10, 
the lowest value in our sample. 

We carry out a similar analysis for Zl and Z2wt to that 
for the R5 object presented in Figure 1. We identify the 
innermost 100 SPH particles at the final simulation output 
and then trace these particles back to earlier times. The time 



evolution of the median values of various physical quantities 
in the Zl and Z2wt simulations are shown in Figure 2. 

On the whole, the formation paths of protostars in the 
Zl and the Z2wt simulations are similar to those in R5. The 
baryon fraction in the halo is initially lower than the uni- 
versal value but it then catches up, as in R5. However, the 
mass scale where this happens is smaller because the cos- 
mological Jeans mass, Mj(z) oc T 3/2 p~ 1/2 oc (1 + z) 3/2 , is 
smaller at lower redshifts. After this stage, the gas is con- 
tinuously compressed and its temperature raises, resulting 
in the rapid formation of H2. When the gas temperature 
reaches T ~ 1000 K and the H2 fraction increases to a few 
times 10 -4 , the gas begins to cool rapidly, just as in R5. 

Despite the overall similarity of the formation paths of 
the first objects in Zl, Z2wt and R5, the timescale for col- 
lapse is very different. In R5 it takes only about 3 million 
years between the time when efficient cooling starts (sig- 
nalled by the decline in temperature) to the time when 
the number density reaches nu ~ 10 10 cm -3 , but this phase 
takes 10 and 100 million years in Zl and Z2wt respectively. 
The initial collapse of protostellar objects proceeds much 
more rapidly at very high redshift because of the higher den- 
sity which accelerates the formation of molecular hydrogen, 
enhances initial cooling, and shortens the dynamical time. 

These results suggest a relatively simple picture of the 
formation path of the first stars although the detailed evo- 
lution can be complicated and may differ from object to 
object depending on the particular infall pattern. As the 
halo virial temperature approaches a critical temperature, 
T ~ 1000 K, the rapid production of molecular hydrogen 
starts once the H2 fraction reaches a critical value of a few 
times 10 -4 , allowing the gas to cool efficiently. The timescale 
for the subsequent evolution to a protostar varies systemat- 
ically depending on the redshift when the collpase is taking 
place, with halos at higher redshift collapsing faster. In our 
simulations the timescale varies by a factor as large as 20 
between redshifts z ~ 50 and z ~ 10. 

Panel (a) in Figures 1 and 2 also suggests that the virial 
temperature provides only a crude estimate of the gas tem- 
perature in the halos because the gas accretion process in 
the protostellar clouds is rather complicated. The maximum 
temperature, which marks the site of the accretion shock, 
can often differ significantly from the virial temperature. 
These two temperatures only roughly match at the time 
when the cosmic Jeans mass is reached, which is also the 
time when the baryon fraction in the halo becomes close to 
the universal value. After this, the maximum temperature 
is usually smaller than the estimated virial temperature by 
a factor as large as 50%. 

Our results imply that a criterion for deciding whether 
or not a halo hosts a star-forming cloud based purely on the 
mass or virial temperature of the host halo, as used exten- 
sively in the literature, is unlikely to give the correct answer. 
While it is true that the gas starts to collapse and cool when 
the halo attains a virial temperature close to ~ 1000 K, the 
processes eventually leading to the formation of a protostar 
can take a significant amount of time and this must be taken 
into account. 
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Figure 2. Evolution of various properties of the final innermost 100 SPH particles in the Zl simulation (left panels) and the Z2wt 
simulation (right panels). Panel (a): the black solid lines with symbols show the virial mass of the main progenitor; the (red) solid lines 
show the mass of the gas component within the virial radius of the halo; the (red) the (red) dashed lines are the fiducial gas mass 
that the halo would contain at the universal baryon fraction, fa = fl^/flm. Panel (b): the evolution of the H2 fraction of the final 100 
innermost gas particles. Panel (c): the gas temperature evolution. The lower (blue) solid line is an estimate of the virial temperature of 
the main halo, as described in the text; the black solid lines are the maximum temperature defined as the maximum of the spherically 
mass-averaged temperature profile of the the gas in the main halo. The (red) solid lines with symbols show the evolution of the median 
value of the temperature for the 100 selected particles. Panel (d): the solid lines with symbols are the median values of the number 
density of hydrogen nuclei in the 100 selected particles. The dashed lines indicate 200 times the cosmological mean number density. 



4.3 Star formation in the R5 simulation 

4-3.1 The distributions of dark matter and gas 

The distributions of dark matter density, gas density and 
mass-weighted temperature for the R5 object are shown, 



at three different epochs, in Figure 3. The images are of 
length 4r2oo on a side and are projected over a depth of 
r2oo centred on the potential minimum. Here r2oo is defined 
as the virial radius of the dominant object at the final time. 
The scales on the axes give the length in physical units and 
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Figure 3. Projected dark matter density (left), gas density (middle) and gas temperature (right), weighted by mass at three redshifts 
in the R5 simulation. In all cases, the projected regions have dimensions of 4 times the virial radius of the R5 object in the plane of the 
plot and are projected over a depth of 1 virial radius at the final time. 



a consistent colour table has been used within each column. 
The three epochs shown are representative of the system 
at the time when its mass is below the Jeans mass, close 
to the Jeans mass, and at the end, when runaway cooling 
is occurring at the centre. From left to right, the columns 
give the dark matter surface density in physical units, the 
gas density in physical units, and the mass- weighted gas 
temperature, respectively. 

At redshift z = 60, the dark matter structure is very 
filamentary, as is typical in CDM models (e.g. Abel, Bryan 
& Norman 2000; Yoshida et al. 2003; Gao et al. 2005). How- 
ever, the gas is noticeably less clumpy than the dark matter. 
Nonetheless, the gas does get compressed due to the signifi- 
cant concentration of the dark matter. As a result, the entire 
region has a temperature of ~ 300 K, several times the tem- 



perature of rsj 80 K expected for the cosmic mean density at 
that epoch. 

At redshift z = 52, the dark matter structures are al- 
ready much more pronounced than at z = 60. Following the 
rapid growth of the dark matter structures, the gas begins to 
settle into the dark matter halos; two dense gas clumps are 
clearly seen in the density map. The increase in gas density 
results in substantial heating of the gas due to compression. 
This can be seen in the large region around the filaments 
where the temperature is ~ 600 K. In the central region of 
the main dark matter halo, two small brighter regions are 
visible with temperature of ~ 1000 K. 

The bottom panels show the system at the final time, 
when star formation is beginning in the innermost regions. 
Note that a second concentration is also evident at a distance 
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of about 40 proper parsecs in projection. In the bottom-right 
panel, two shocks develop along the filaments. On a scale 
of ~ 10 pc, the shock around the dominant object is even 
more pronounced. Note that the shock occurs in the very 
inner regions rather than at the viral radius, r v i r = 35 pc. 
Just inside the shock front, a molecular hyrdogen rich and 
rapidly cooling region is found. 

A more detailed set of images of the R5 object at the 
final time is displayed in Figure 4. From left to right, the 
columns show surface density, mass- weighted temperature, 
and the radial velocity with respect to the star-forming core. 
From top to bottom, the images show successive zooms by 
factors of 8 in linear scale onto the centre of the gravitational 
potential at the final time. The radial velocity is defined as 
the velocity relative to the average velocity of the innermost 
100 SPH particles. On small scales, the gas cloud reveals a 
complex structure, especially in the radial velocity map. The 
gas accretion is clearly asymmetrical. The infall velocity is 
much larger for the gas on the right hand side, where the 
secondary gas concentration is closing in towards the dom- 
inant object, accompanied by a one-sided accretion shock 
at a scale of ~ 1 pc. Clearly, describing these processes as 
spherically symmetric gas accretion is unrealistic. 

4-3.2 Time evolution of radial profiles around the 
dominant object 

The time evolution of the gas in the R5 object is described 
in Figure 5 which shows the mass- weighted radial profiles of 
various physical quantities of interest. We define the cool- 
ing timescale as t coo i = — dt/dlnT, where T is the gas 
temperature, and t is time. The free-fall time is defined as 
tff = 27ry3/32Gp, where p is the density. We show 6 epochs 
corresponding to z = 52.24, 50.29, 48.50, 48.02 and 47.91; 
the final epoch corresponds to z = 47.90. 

We start 6.4 million years before the final time, at 
z = 52.24, when the R5 halo has only just reached the Jeans 
mass. Gravity is then strong enough to cause gas to set- 
tle into the potential well of the halo. By then, the central 
part of the gas distribution is heated to more than 1000 K. 
However, the cooling time is still much longer than the free- 
fall time because the gas density and molecular hydrogen 
fraction are still relatively low. At this stage, the gas con- 
tracts slowly. About 2.7 million years later, at z = 50.29, 
the central gas density has increased by a factor of 10, and 
the molecular fraction has also increased by a factor of 10. 
This results in a much shorter cooling time for the central 
2000 M© of gas, which is now comparable to the free-fall 
time. At this point, the gas cloud begins to collapse rapidly. 

Another 2.7 million years later, at z — 48.50, the central 
1000 Mq of gas approach the critical density, nu ~ 10 4 cm -3 , 
which marks the transition of H2-cooling from NLTE rota- 
tional level populations to an LTE state. Above this den- 
sity, the cooling rate per molecule of hydrogen becomes in- 
dependent of density. The central gas temperature falls to 
close to 200 K, which is roughly the minimum temperature 
that can be reached by molecular hydrogen cooling. By now 
the production of molecular hydrogen has entered a regime 
where the abundance scales logarithmically with time due to 
the continuous depletion of electrons through recombination 
(Tegmark et al. 1997). 

Finally, a further 1 million years later, at z = 48.02, the 



central gas density approaches the critical density required 
for three-body reactions, resulting in rapid H2 formation and 
an increase of the gas temperature due to heating from both 
compression and the release of molecular binding energy. 
After 10 5 years, there is a substantial increase in the H2 
fraction in the central few hundred solar masses of gas with 
density exceeding 10 8 cm -3 . As a consequence, the cooling 
time of the central 20 M of gas becomes shorter than the 
free-fall time, and the mass eventually exceeds the locally 
estimated Bonnor-Elbert mass (Elbert 1955; Bonnor 1956): 

M BE ~ 2OM T 3/2 n" 1/2 /i" 2 7 2 . (3) 

Runaway collapse is triggered at this time. Here n is the 
number density of hydrogen nuclei, p is the molecular weight 
and 7 is the adiabatic index. Note that the Bonnor-Elbert 
mass is similar to the Jeans mass but it assumes an isother- 
mal rather than an uniform gas distribution. 



5 OBJECT-TO-OBJECT COMPARISONS: 
DIFFERENCENS AND SIMILARITIES 

As we have seen in the previous section, the formation epoch 
of the first stars can span a broad redshift range, and there 
are systematic differences in the formation process of proto- 
stellar clouds that depend on redshift. An additional source 
of variation comes from the differences in the assembly his- 
tories of these objects at a given redshift (e.g. Lacey & Cole 
1993). It is therefore important to analyse a number of real- 
isations in order to assess the expected scatter in the prop- 
erties of star-forming clouds. We now address this issue by 
jointly considering all of our simulations. 

5.1 Radial profiles 

Figure 6 shows the radial profiles of the density, enclosed 
gas mass, molecular hydrogen fraction and temperature in 
all our simulations. The comparison is made at the final 
time when the central gas has roughly reached the limiting 
density, nu ~ 10 10 cm -3 , in all the simulations. At this time, 
runaway collapse has set in and the minimum simulation 
timestep becomes extremely small. The top left panel shows 
the number density profile of hydrogen nuclei. Each object 
is represented by a different colour line as indicated in the 
top-left panel. A striking feature in this plot is that the gas 
density profiles of the objects in all our simulations coincide 
with one another on scales less than ~ 10 pc, in spite of their 
very different collapse redshift s and environments. 

The scatter is remarkably small, only a factor of ~ 2 
about the mean, which is smaller than in the simulations of 
O'Shea &Norman (2006b). For comparison, the red crosses 
and the blue dashed line show the fiducial density profiles of 
ABN02 and Yoshida et al. (2006), respectively. Our simula- 
tions agree well with ABN02 on scales smaller than 0.02 pc, 
but not very well at larger radii. As we will show in Sec- 
tion 6, the different baryon fractions assumed in the two 
studies largely account for this discrepancy. Yoshida et al. 
(2006), on the other hand, adopt the same cosmology as us. 
Our results agree well within 0.2 pc, but not as well on larger 
scales. This suggests that the actual variation in density pro- 
file from object to object may be bigger than indicated by 
the limited samples we have simulated here. The top right 
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Figure 5. Radially averaged profiles of various physical quantities at different epochs in the R5 simulation. The colour lines refer to the 
epochs indicated in panel (b), where the time is expressed as the number of years before the final time, and the corresponding redshifts 
from top to bottom are: 47.90 (the final time), 47.91, z = 48.02, z = 48.50, z = 50.29, and z = 52.24. Panel (a): enclosed gas mass as 
a function of radius. Panel (b): number density profile of hydrogen nuclei. The dashed lines are the local Bonnor-Ebert mass, which is 
similar to the Jeans mass but assumes an isothermal rather than an adiabatic distribution for the gas. Panel (c): molecular hydrogen 
fraction as a function of the enclosed gas mass. Panel (d): gas temperature. Panel (e): the ratio of free-fall time to the cooling timescale. 
The vertical doted lines in all panels indicate the virial radius of the halo at the final time. 



panel shows the enclosed gas mass as a function of radius, 
which also confirms the rather small variation in mass profile 
from object to object in our simulations. 

The molecular hydrogen fractions in all our simulations 
also agree well with each other at radii less than 10 pc. The 
temperature structures are similar as well. The latter can 
be roughly divided into 4 physically distinct zones. Zone 
1 is the outermost region and corresponds to cosmic in- 
fall, where the gas is continually heating up, culminating at 
the accretion shock front where the maximum temperature 
occurs. In the region immediately behind the shock front, 
which makes up zone 2, the gas temperature decreases in- 
wards up to a radius, ~ 1 pc, where the temperature reaches 



a minimum of ~ 200 K, below which the molecular hydro- 
gen cooling rate drops exponentially. Within this radius, in 
zone 3, there is a second infall region in which the tem- 
perature rises back up again and approaches 1000 K at a 
radius of ~ 0.01 pc. Here, the gas density reaches a value 
of 71h ~ 10 8 cm -3 . At this characteristic density, three-body 
reactions which rapidly form molecular hydrogen start to 
operate and the enhanced cooling associated with the rising 
molecular hydrogen fraction causes the temperature to de- 
crease towards smaller radii temporarily. The temperature 
profiles are all qualitatively similar but they can be grouped 
into two classes: R5, R5wt, Zlw and Z2wt in one, and Zl, 
Zlwt, Z2, Z2w in the other. As we will see in the next sub- 
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Figure 6. Radially averaged profiles of various physical quantities as a function of radius in the 8 simulations. The different lines refer 
to the different simulations as indicated by the labels in the top left-hand panel. Top left: number density of hydrogen nuclei. Top right: 
enclosed gas mass. Bottom left: molecular hydrogen fraction. Bottom right: gas temperature. 



section, the star- forming clouds in the second group have a 
disklike structure. As a result, the collapse is slower and the 
competition between compressional heating and molecular 
hydrogen cooling results in a lower inner temperature. 

Figure 7 compares the radial profiles of various dynam- 
ical quantities. Here, for clarity, we show results only for 3 
typical simulations, R5, Zl and Z2wt, which we analysed in 
Section 3. A complete version of the figure including all the 
other simulations may be found in the Appendix. Figure 7 
reveals that despite the similarities in the gas properties seen 
in Figure 6, the internal gas motions differ dramatically from 



object to object. They are only similar in the cosmic infall 
region which terminates at the first accretion shock where 
the radial velocity (defined to be negative for infalling gas) 
reaches a local minimum and the temperature a local max- 
imum. The mass inside this radius is about 10 4 Mq for the 
R5 and Zl objects, and 10 5 M for the Z2wt object. For the 
R5 object, the gas infall velocity is relatively constant but 
declines rapidly at a radius corresponding to a mass scale 
of 1O 3 M . Just inside this scale, the gas infall speeds up 
again. For the Zl object, the infall speed increases just be- 
hind the accretion shock and reaches a local maximum at 
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Figure 7. Radial profiles of various physical quantities as a function of the enclosed gas mass at the final time in 3 representative 
simulations: R5, Zl and Z2wt. Panel (a): profiles of negative tangential velocity. Panel (b): radial velocity profiles. Panel (c): the ratio 
of the free-fall time to the cooling time. Panel (d): the ratio of the tangential velocity to the Keplerian velocity v Kepler — (GM/r) 1 / 2 . 
Panel (e): the orientation of the angular momentum. Note that only the angle <p is shown. 



a mass scale of 3 x 10 3 M©, where the ratio of the free-fall 
time to the cooling time is maximal (panel c) . Exactly at the 
same scale, the tangential velocity starts to become larger 
than the infall velocity. For the Z2wt object, the radial ve- 
locity stays constant at large radii. At a mass scale of about 
5 x 10 3 M©, the radial infall velocity increases until a mass 
scale of a few tens of solar masses. Note that a jump in the 
ratio of the free-fall time to the cooling time often marks 
a transition in the velocity and orientation of the angular 
momentum. 



Panel (d) gives the ratio of the actual tangential velocity 
of the gas to the velocity required for Keplerian motion, de- 
fined as |14|/Vkepierian- The curves for the three simulations 
differ substantially, implying that the degree of rotational 
support is rather different. Finally, panel (e) shows the az- 
imuthal angle between the z-axis in the simulation and the 
angular momentum vector of the gas interior to a given ra- 
dius, or equivalent ly an enclosed mass. It is clear that the 
different parts of the halos often spin in different directions. 
For the R5 object, the gas only rotates in the same direction 
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inside an enclosed mass of 300 M©, while the same is true 
for the Zl and Z2wt halos inside 3000 M and 10000 M , re- 
spectively. Comparing with the other panels in the Figure, 
we see that a jump in tangential velocity typically accom- 
panies a jump in the orientation of the angular momentum. 

5.2 Morphology of the star- forming clouds 

In spite of the similarities in radial profiles seen in the previ- 
ous subsection, the final star- forming clouds exhibit a vari- 
ety of morphologies. We define the edge of the star-forming 
cloud as the largest radius inside which the orientation of the 
angular momentum of the enclosed material (determined by 
the azimuthal angle introduced in the previous subsection) 
does not change significantly. We then calculate the prin- 
cipal axratios, q — b/a and s — c/a (taking a > b > c), 
of the inertia tensor of all gas particles in the star-forming 
cloud. (These axial ratios may be viewed as the equivalent 
axial ratios of a homogeneous ellipsoidal distribution with 
the same inertia tensor as the gas distribution.) The axial 
ratios for all our simulated objects are given in Table 5. We 
also list here the ratio of the tangential velocity to the Ke- 
plerian velocity at the edge of the cloud. It is clear that the 
star-forming clouds are rarely close to spherical. Only three 
objects in our sample, R5, R5wt and Zl, have a relatively 
round shape, with a tangential velocity that is only ~ 30% of 
that required for Keplerian motion. These clouds are similar 
in shape to the objects simulated by ABN02 and Yoshida et 
al. (2006). However, half of our objects, Zl, Zlwt, Z2 and 
Z2w, have a disk-like structure of varying size and thickness 
which is nearly rotational supported {Vt /Vxepierian > 0.6), 
while one, Z2Wt, has a bar- like structure. 

Density and temperature maps at the final time in the 
Zl and the Z2wt objects are shown in Figures 8 and 9 re- 
spectively. The left and right panels give the density and 
temperature maps respectively. The top panels depict square 
regions of length 4 x r2oo and depth r2oo in projection. In 
the lower panels, the scale has been adjusted to illustrate the 
disk in Zl and the bar in Z2wt, with the projected velocity 
field superposed onto the density map. Both structures are 
shown face- on. 

There are systematic differences between the internal 
structure of the main object and its surroundings in R5, Zl 
and Z2wt. For the R5 object, the surrounding dark matter 
distribution is quite complex. There are 3 prominent objects 
connected to each other by filaments (see Fig 3) and the 
gas at the intersection of these filaments is shock-heated 
by infall. For the Zl and Z2wt objects, the arrangement 
of the dark matter is simpler. These objects form at the 
intersection of only two filaments. For the gas distribution, 
the situation is the reverse: the gas around R5 is relatively 
featureless, but in Zl and Z2wt, the inflowing gas has a 
complex structure. This is reflected also in the temperature 
maps of the Zl and Z2wt objects which show a distinctly 
chaotic pattern. 

The star-forming cloud in the R5 object is relatively 
spherical, as may be seen in Table 5, while the Zl object has 
a disk- like structure extending out to a radius of ~ 0.1 pc. 
The Z2wt object has a huge bar extending out to a radius 
of rsj 1 pc. In both these objects the gas velocity fields are 
well organized, revealing a nearly circular pattern in Zl and 
a vortical pattern around the bar in Z2wt. However, in the 



very inner regions of Zl, the gas has less rotational support 
(see Fig. 7) and is undergoing an inside-out collapse. It will 
be very interesting to study the fate of this gas in future 
work. 

Disk-like structures have been seen in previous simu- 
lations of first stars by Bromm et al. (1999, 2002). These 
authors used constrained initial conditions and arbitrarily 
assigned a high value of the spin parameter to the simulated 
region. The disk in our simulation is substantially differ- 
ent from theirs. In Bromm et al. (1999, 2000), the disk- like 
structure formed shortly after virialization, then broke up 
rapidly into several disjointed gas clumps ranging in mass 
from 200 to 1O 4 M . In our simulation which started from 
proper cosmological initial conditions, the gas acquires an- 
gular momentum self-consistently by tidal torques from sur- 
rounding structures, rather than being put in by hand as in 
Bromm et al. The disk forms at a later stage when the cen- 
tral gas density approaches 10 9 cm -3 , and there is no sign 
that the disk- like structure is unstable to fragmentation up 
to the time when we stop our simulation. This is because 
the disk is hot and thick, and so is pressure supported. We 
have checked that the Toomre Q parameter, 



of the Zl disk lies in the range 2 — 3 which also suggests 
that it is stable against perturbations. Here, Q is the local 
angular frequency, c s is the sound speed, and E is the surface 
density. 

The Zlwt and Z2w simulations also formed thick disks 
and careful inspection reveals the presence of spiral arms, as 
may be seen in Fig. 10. The associated asymmetric density 
fluctuations will produce gravitational torques that will help 
transport angular momentum outwards efficiently. 

The processes that drive the different morphologies ex- 
hibited by the star-forming clouds in our simulations remain 
unclear. We have been unable to find any correlation be- 
tween morphology and the various global properties that 
we have examined and we have found no clear relationship 
to the merging history of the host dark halo. We defer a 
more detailed study of the origin of the morphologies of 
star-forming clouds to a future study. 

5.3 Accretion rate 

The rate of mass accretion within the pre-stellar cloud is 
an important quantity which governs the evolution of the 
protostar. The growth of the protostar can be estimated by 
the mass accretion rate, M, via 

M*(t) =M prot o+ f M(t)dt, (5) 
Jo 

where M pro to is an initial protostellar mass. In our simula- 
tions, we have neither enough resolution nor a sufficiently 
complete physical model to be able to resolve the initial 
protostar. One-dimensional hydrodynamic simulations in- 
cluding all the relevant physics suggest that the protostar 
mass is of order a few Jupiter masses (Omukai& Nishi 1998; 
Ripamonti et al. 2002). 

Ideally, we need to follow the protostellar evolution 
along with the accretion in order to obtain a reliable fi- 
nal mass for the star. However, the coupled evolution of the 
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s = c/a 
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Table 5. The principal axes ratios, q = b/a and s = c/a, for our 8 simulations, where a > b > c. Here, the edge of the star-forming cloud 
is defined as the outer radius of the region in which the orientation of the specific angular momentum is roughly constant. This radius 
is denoted by R CO re] M core is the enclosed gas mass within it. 




-150 -100 -50 50 100 1 50 -150 -100 -50 50 1 00 150 

x[pc] x[pc] 




-0.2 -0.1 0.0 0.1 0.2 -0.2 -0.1 0.0 0.1 0.2 

x[pc] x[pc] 



Figure 8. Density (left) and temperature (right) maps for the zl simulation which has a disk-like structure at the final time. The upper 
panel shows a slice of side Ar v i r and depth r v i r . The bottom panel shows a cubic region of side 0.4 pc chosen to include the disk-like 
structure which is shown face-on. The projected velocity field is superposed onto the density map. 
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Figure 9. Density (left) and temperature (right) maps for the Z2wt simulation which has a bar-like structure at the final time. The 
upper panel shows a slice of side Ar v i r and depth r vir . The bottom panel shows a cubic region chosen to include the bar-like structure 
which is shown face-on. The projected velocity field is superposed onto the density map. 



protostar and the accreting envelope gas is a very complex 
physical problem which has not been solved so far, even 
in one-dimensional radiative hydro-dynamical simulations. 
At present, to estimate plausible stellar masses one needs 
to resort to either analytical collapse solutions (e.g. Larsen 
1968; Penston 1967, Shu 1977), to detailed protostellar cal- 
culations with constant accretion rates (Omukai & Palla 
2001;2003), or to estimates based on the instantaneous ac- 
cretion rate from a particular numerical simulation (Tan & 
McKee 2003; Omukai & Palla 2003; Yoshida et al. 2006; but 
see Bromm & Loeb 2004 for a slightly different approach). 
Here we present the instantaneous accretion rates found in 
our own simulations. 

In Figure 11, we show the instantaneous accretion rate, 
as a function of enclosed gas mass, for our 8 simulations. 



Here, the accretion rate is defined as M(r) = Aitr 2 p(r)\v Ta d\ 
and is calculated at the time when the central gas density 
has reached a value of nn ~ 10 10 cm -3 . Since the accretion 
rate depends on the local gas density and on the gas infall 
velocity, we expect a considerable scatter from one object to 
another due to very different gas motions and variations in 
the local density. 

The R5, R5wt, Zlw and Z2wt objects have the largest 
instantaneous accretion rates amongst all the simulations, 
and, at the same time, they also have the least rotational 
support. The Zl, Zlwt, Z2, and Z2W objects have the small- 
est accretion rates because their disk-like structures slow 
down the gas infall rate. Thus, the accretion time allows 
our simulated objects to be roughly divided into two pop- 
ulations, depending on the morphology of the star-forming 
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Figure 10. Density maps for the Zlwt and the Z2w simulations which have a disk-like structure at the final time. These images clearly 
reveal the presence of spiral arms contained in the disk-like structure. The projection is in an arbitrary direction. 
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Figure 11. Accretion rate, M(r) = Attt 2 p(r)v[ n f a i\, as a function 
of enclosed gas mass for our 8 simulations. Here, the accretion rate 
is computed at the time when the central gas density has reached 
a value of tih ~ 10 10 cm -3 . The accretion rates for objects with a 
central disk-like structure are ploted as dahsed lines, while objects 
with other morphologies are plotted as solid lines. 



clouds. The accretion time for rotationally supported disks 
is much longer than for the other morphologies. 

Omukai & Palla (2003) identified a critical accretion 
rate of ~ 4 x 10 -3 M©yr _1 (the thick solid line in Figure 11) 
for a star to form. This rate is physically motivated by the 
fact that the protostar luminosity cannot exceed the Ed- 
dington limit; otherwise a radiation-driven expansion of the 
protostar would halt or even reverse the accretion when the 



star reaches a mass of ~ 100 M©. Note that with this high 
accretion rate, the star may in fact start burning hydrogen 
well before the accretion slows down, so that M/M becomes 
much longer than the lifetime of the star. Such an object 
would not be a zero-age main sequence mass. If the accre- 
tion rate is smaller than the critical rate, however, the star 
could continue to accrete gas after it settles onto the zero- 
age main sequence (ZAMS), allowing it to reach a larger final 
mass. In our simulations, objects without disks all have in- 
stantaneous accretion rates on scales corresponding to an 
enclosed mass of less than a few hundred solar masses that 
are a few times higher than the critical rate, while those 
with disk-like star-forming clouds have lower or comparable 
rates to the critical one. 

Does this mean that our simulated objects with accre- 
tion rates larger than the critical value would be limited to 
a final mass of about ^100 M© by radiative feedback from 
their own protostars? To address this question, we carry 
out an experimental calculation to examine whether the in- 
stantaneous accretion rate is robust enough to serve as a 
reliable estimator of the theoretically expected final mass 
of the star. To this end, we examine whether the velocity 
structure of the gas in the protostellar environment changes 
rapidly, particularly whether the gas would spin up after 
further collapse which may slow down (or even halt) later 
accretion. For the R5 simulation, which has a star- forming 
cloud with little rotational support at our fiducial end state, 
we continued to integrate forward in time for another 1000 
years (unlike ABN02 and Yoshida et al. (2006) who halted 
their simulation at this time). The ratio of the tangential 
velocity to the Keplerian rotational velocity at several times 
during this period is shown in Figure 12. 

The red solid line with open diamonds refers to the final 
time when our physical model is still adequate. One can see 
that as the evolution is allowed to proceed, the gas cloud be- 
comes increasingly rotationally supported and forms a ring- 
like structure at a length scale of ~ 0.02pc within one thou- 



© 0000 RAS, MNRAS 000, 000-000 



18 Gao L. et al 



sand years, which is a short time in the lifetime of a proto- 
star. While this evolution is not fully realistic because nei- 
ther is the adopted inner boundary condition correct nor is 
the feedback generated by the protostar considered, it never- 
theless highlights the fact that the velocity field in the proto- 
stellar environment can change rather rapidly. This suggests 
that it is hard to draw definitive conclusions regarding the 
likely mass of the star-forming cloud in the R5 halo based 
solely on the instantaneous accretion rate. 

Continuing the integration of the other systems which 
have similar morphologies to R5 yields similar results: the 
accretion rate typically tends to become smaller as gas be- 
comes more and more rotationally supported. By contrast, 
in the Z2wt object, whose star- forming cloud has a bar- 
like structure, the accretion rate rises at slightly later times 
when the asymmetrical mode of the bar transfers angular 
momentum outwards very rapidly and so boosts gas infall. 
As a result, gas accretion onto the protostar of the Z2w 
halo might be limited by the feedback effects discussed in 
Omukai & Palla (2003) and Tan & McKee (2003). Integrat- 
ing the systems with a disk-like structure further in time, we 
find that the accretion rate falls even below that measured 
at the final reliable output time. We conclude that the in- 
stantaneous accretion is not a robust indicator of the mass 
of the final PoP-III star. 

The instantaneous accretion rates within a radius con- 
taining 100 solar masses found in the simulations by ABN02 
and Yoshida et al. (2006) are compatible with our results. 
Their rates lie within the scatter we find from our set of 
eight simulations. Recently, O'Shea & Norman (2006b) car- 
ried out a set of simulations similar to ours although using 
an AMR code. While they found, as we do, that there is a 
considerable scatter in the instantaneous accretion rate from 
object to object, their scatter is an order of magnitude larger 
than ours. Another significant difference between their work 
and ours is that we find no correlation between the redshift 
of formation and the instantaneous accretion rate for our ob- 
jects. We do not know the reasons for these differences and 
it will be interesting to run the simulations with identical 
initial conditions with both our code and Enzo. 



6 INFLUENCE OF THE BARYON FRACTION 

It is interesting to examine how our choice of cosmological 
parameters affects the properties of primordial gas clouds. 
In particular, the baryonic matter density, or baryon frac- 
tion, may be one parameter which can signficantly affect the 
gas evolution. Naively, one would expect that the baryon 
fraction influences the gas collapse process, since both the 
production of molecular hydrogen and the cooling function 
at low densities depend strongly on density. Furthermore, 
initially the gas distribution in dark matter halos follows 
closely that of the dark matter component, with an enclosed 
baryonic fraction that should be close to the mean cosmic 
value (e.g. Lin et al. 2006). Thus, varying the baryon frac- 
tion may significantly change the density distribution of the 
gas within a halo. 

Given the uncertainties that remain regarding the pre- 
cise values of the basic cosmological parameters of our Uni- 
verse (see the discussion in Spergel et al. 2006), in this sec- 
tion, we investigate the effects of varying the baryon fraction, 
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Figure 12. The ratio of the tangential velocity to the Keplerian 
circular velocity, as a function of radius, at several epochs for the 
R5 simulation. The red curve with open squares refers to our last 
trusted epoch when the optically thin assumption is still valid. 
The time intervals between the other outputs and this final time 
are indicated in the labels. 



f b = Q b /Q m , on the formation of the first pro-stellar objects. 
The precise values of the matter-density, Q m , of the dark en- 
ergy component, Q\, and the shape and normalization of the 
initial power spectrum are all important for determining the 
abundance and the formation epoch of the first stars, but 
are less important than the baryon fraction in determining 
the properties of individual star forming clouds. We consider 
some of these other properties in the next section. A final 
reason for varying the baryon fraction is that it facilitates 
a direct comparison with the results of ABN02 who used a 
smaller value than we have assumed so far. 

We have resimulated the Zl object with three different 
choices for the baryon fraction. Our fiducial choice is fb = 
2/15, until recently, the value assumed in the standard (or 
"concordance") ACDM model. We compare this with the 
lower baryon fraction, fb = 0.06, which was often used in 
the context of the now abandoned Q m = 1 standard CDM 
model assumed by ABN02. Finally, we also consider fb = 
0.178, as suggested by the latest WMAP and 2dFGRS results 
(Sanchez et al. 2005, Spergel et al. 2006). We refer to the 
simulations with lower and higher baryon fractions as as 
Zllb and Zlhv respectively. 

In Figure 13, we compare the radial profiles of the gas 
density, enclosed gas mass, temperature and radial velocity 
in the three simulations. The magenta dashed lines show 
the results of ABN02. The radially averaged physical prop- 
erties of the Zlhb and Zl simulations are generally very 
similar, although some small differences are apparent. For 
example, the density profile of Zlhb is ~ 20% higher than 
that of Zl over much of the range shown. The differences 
are much larger for the lower baryon fraction simulation. 
On scales 0.1 — 1 pc, the gas density in Zllb is roughly a 
factor of 2 lower than in the other simulations. The radial 
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Figure 13. Comparison of various spherically averaged profiles in the Zl simulation for different assumed baryon fractions. Panel (a): 
enclosed gas mass; panel (b): number density of hydrogen nuclei; panel (c): temperature; panel (d): radial velocity. The lines refer to 
different simulations as indicated by the label. The redshifts are z = 25.93 for the Zl run, z = 26.26 for the Zlhb simulation, and 
z = 24.26 for the Zllb run. Relative to the fiducial, fa = 2/15, case, star formation in the higher baryon fraction simulation occurs ~ 2 
million years earlier, while in the lower baryon fraction simulation, it is delayed by ~ 12 million years. 



gas velocity of Zllb is also very different from that of Zl and 
Zllb. In fact, surprisingly, the disk- like structure has disap- 
peared in the Zllb simulation, although it is still present 
in the higher baryon fraction simulation, Zlhb. Finally, we 
note that matching the lower baryon fraction assumed by 
ABN02, as in Zllb, produces much better agreement be- 
tween our results and theirs. 

We find that the baryon fraction influences the collapse 
time of the star- forming cloud quite significantly. Relative to 



Zl, the collapse proceeds slightly faster in Zlhb and much 
more slowly in Zllb. In Zlhb, star formation occurs as early 
as z — 26.26 (that is 2 million years earlier than in Zl), 
while in Zllb, it is delayed by 12 million years, to a redshift 
of z = 24.26. We conclude that a slight change in the baryon 
fraction does not affect the properties of star-forming clouds 
much, but changes by a factor ~ 2 have a significant effect. 
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Figure 14. Left-hand panel: halo mass as a function of redshift. The solid lines corresond to the mass M*, denned through u = cr(M)/5 c = 
1; the dashed lines correspond to the mass of 2cr halos {y = 2); the dot-dashed lines correspond to the mass of 3cr (y = 3) halos. The 
black curves show results for cosmological parameters derived from the 2dFGRS and the first year WMAP data and the red curves for 
parameters derived from the 2dFGRS and three-year WMAP data (Sanchez et al. 2005, Spergel et al 2006). The two long-dashed blue 
lines show the halo mass which, at each redshift, has a virial temperature of more than 1000 K and 10000 K, respectively, as indicated by 
the labels. Right-hand panel: the predicted halo abundance in the two cosmologies. The solid lines show cumulative comoving number 
densities of halos more massive than a given mass as a function of redshift (upper axis). The dashed lines show the predicted comoving 
number density of halos that have a virial temperature exceeding 1000 K as a function of redshift (lower axis). This is roughly the critical 
temperature at which gas can start efficient cooling by molecular hydrogen. 



7 ABUNDANCE OF STAR-FORMING HALOS 

So far, we have concentrated on the physical processes that 
lead to the formation of the first protostars in a ACDM 
universe. We now turn to a discussion of the abundance of 
early star-forming halos and its dependance on cosmologi- 
cal parameters in a ACDM universe. As we have shown in 
Section 3, the basic principles that determine the formation 
paths of the first protostars are simple: once gas is heated 
up to a temperature of around 1000 K, the production of 
molecular hydrogen is boosted sufficiently for the gas cloud 
to cool very rapidly. The ensuing collapse proceeds faster 
if it takes place in a halo at higher redshift because in that 
case the gas starts off with a higher density and temperature 
and a shorter dynamical time. 

As shown in Figures 1 and 5, the virial temperature 
provides a reasonable estimate of the initial temperature of 
the star-forming material at the time when it begins to settle 
into the halo. In Figure 14, we show the abundance of objects 
with virial temperature above 1000 K. Note that the gas in a 
halo with this temperature is predicted to start collapsing at 
this time, but the collapse process itself, eventually yielding 
a protostellar object, takes a further ~ 10 6 — 10 8 years. This 
corresponds to is significant interval in redshift, Az ^ of a 
few, depending on the exact redshift. 

The first stars in the standard ACDM cosmology are 
often assumed to form typically at redshifts between z = 20 
and z = 30 because 3a halos with a virial temperature of 
c± 1000 K collapse in this redshift range. In the left-hand 



panel of Figure 14, we plot the abundance of la, 2a and 
3cr halos as a function of redshift. A similar figure has been 
shown in a number of previous studies (e.g. Barkana & Loeb 
2001) but here we compare results obtained for two choices 
of cosmological parameters, those inferred by combining the 
galaxy power spectrum of the 2dFGRS with the first-year 
and three-year WMAP data respectively. In the figure, la 
halos are shown with solid lines, 2a halos with short-dashed 
lines and 3a with dot-dash lines. The two long-dashed blue 
lines show the halo mass which, at each redshift, has a virial 
temperature of more than 1000 K and 10000 K, respectively. 

As a consequence of the significant suppression of power 
on small scales in the 2dFGRS+three-year WMAP cosmol- 
ogy (due to a lower value of as and a red tilt in the pri- 
mordial power spectrum index), the curve for the 3a halos 
in this cosmology is shifted so much that it nearly coincides 
with that for the 2a halos in the older cosmology. The curve 
for halos which have a virial temperature exceeding 1000 K 
crosses the curve for 3a halos in the 2dFGRS+three-year 
WMAP cosmology at z ~ 15 instead of at z ~ 25 as in the 
older cosmology. 

The dramatic changes in the number of halos that can 
host the first stars arising from the relatively small change 
in the values of the cosmological parameters in the 2 models 
we considering are illustrated in the right-hand panel of Fig- 
ure 14. Here, the comoving abundance of potential first-star- 
bearing halos, i.e. those with a virial temperature exceeding 
1000 K is plotted as a function of redshift (labelled along 
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the lower axis) for the two sets of cosmological parameters 
(distinguished by the two colours as in the left panel). For 
reference, we also plot the comoving abundance of all ha- 
los at the present day, as a function of halo mass (labelled 
along the upper axis), based on the fitting formula of Sheth, 
Mo & Tormen (2000), which works reasonably well for large 
masses (Jenkins et al. 2001; Springel et al. 2005; Reed et al. 
2007). 

The abundance of potentially star-forming halos differs 
dramatically in the two cosmologies at redshifts greater than 
z = 10. At redshift z = 50, the predicted abundance of ha- 
los capable of molecular cooling is ~ 7 orders of magnitudes 
lower in the 2dFGRS+three-year WMAP cosmology than 
in the older 2dFGRS+nrst-year WMAP cosmology. As cos- 
mic time increases, the predicted abundance of star-forming 
halos increases very rapidly, as noted by G05. By redshift 
z ~ 35, the abundance of halos with T V - 1T > 1000K has in- 
creased by 6 orders of magnitude, making it comparable in 
the cosmology that uses the 3-year WMAP data to that of 
dwarf galaxies at the present day. By redshift z ~ 15, the 
abundance of such objects is 1000 times higher than that of 
10 8 M© halos today. It is clear that the current uncertainty 
in the values of the cosmological parameters, particularly as, 
is a major source of uncertainty for the predicted abundance 
and evolution of Pop-Ill stars. 



8 SUMMARY AND DISCUSSION 

In this paper, we have presented a large set of cosmologi- 
cal SPH simulations of the formation of the first stars. Our 
simulations reach sub-solar mass resolution and include an 
accurate and detailed chemistry network to treat radiative 
cooling by molecular hydrogen. Our large suite of simula- 
tions allows us to investigate the scatter in the properties 
of star-forming clouds and the redshift dependence of their 
formation process. Our simulations stop when the central 
gas density has reached nu ~ 10 10 cm -3 , the point at which 
our assumption that the gas is optically thin begins to break 
down. 

Despite the complexities in detailed formation paths, 
our simulations suggest a relatively simple picture leading 
up to the development of the first star-forming clouds. Once 
primordial gas is heated to ~ 1000 K, either by shocks or by 
adiabatic compression, the production of molecular hydro- 
gen is boosted. When the molecular fraction reaches a few 
times 10 -4 , the gas cools rapidly and settles to the centre 
of the gravitational potential well of the halo. Our simula- 
tions bear out the same sequence of events sketched in ear- 
lier semi- analytical work by Tegmark et al. (1997) and seen 
in the detailed simulations of Fuller & Couchman (2000), 
Machacek et al. (2001) and Yoshida et al. (2003). Our main 
results may be summarized as follows: 

(1) The formation path of Pop-Ill stars is broadly similar 
in all halos, independently of collapse redshift and environ- 
ment. However, the timescale for the collapse differs sys- 
tematically depending on the redshift at which the collapse 
takes place. The gas condenses more rapidly into stars at 
high redshift and so the formation timescale can differ by 
factors as large as 20 for halos that collapse between red- 
shifts z rsj 50 and z ~ 15. For example, for the R5 object 
which starts to form at z ~ 50, the collpase takes as little as 



a few million years whereas for the Z2wt halo, which began 
to form at z ~ 15, the collapse took 100 million years. 

(2) The spherically averaged radial gas density profiles of 
the star-forming clouds agree well with each other up to a ra- 
dius of rsj 10 pc. The profiles exhibit a power-law distribution 
pocr -2-2-2-3 which extends over 10 orders of magnitude in 
density and is independent of formation redshift and envi- 
ronment. The molecular hydrogen fractions and temperature 
structures are also very similar. The main reason why these 
similarities are present is that once the gas has started to 
collapse, it largely isolates itself from the rest of the evolving 
halo and so the final state of the star-forming cloud becomes 
largely insensitive to the global properties of the host halo 
such as its virial mass, temperature or spin. 

(3) Despite the similarity in the radial profiles of gas den- 
sity, temperature, and molecular hydrogen fraction, the ve- 
locity field of the gas can differ significantly from object to 
object. 

(5) At the moment when our simulations stop, the star- 
forming clouds exhibit a variety of morphologies. Some of 
them are nearly rotationally supported disks, others resem- 
ble flattened spheroids, while others possess self- gravitating 
bars. The instantaneous accretion rate onto the central re- 
gion distinguishes between two kinds of star- forming clouds. 
Those with a disk-like structure have the smaller accretion 
rates, while those that resemble flattened spheroids or pos- 
sess bars have the larger accretion rates. The scatter in ac- 
cretion rate can be more than a factor of 10 and is related 
to the morphology of the cloud. 

(6) A higher baryon fraction increases the gas density ev- 
erywhere within a halo and speeds up the formation of the 
first stars. 

(7) The predicted abundance of Pop-Ill halos in the 
ACDM model varies dramatically depending on whether one 
assumes the cosmological parameters derived using either 
the first year or the three-year WMAP data. For example, 
the abundance of objects similar to our R5 halo differs by 7 
orders of magnitudes at redshift 50 in the two cosmologies. 

At the moment when our simulations stop, i.e. when 
71h reaches a value ~ 10 10 cm -3 , none of the central star- 
forming clouds have fragmented. This is in agreement with 
earlier simulations (ABN02; Bromm et al. 2002; O'Shea & 
Norman 2006a,b; Yoshida et al. 2006). At this time, all our 
simulated objects experience a much higher instantaneous 
accretion rate onto their central region than is common in 
protostars today. If fragmentation does not occur and the 
accretion does not slow down significantly, these objects are 
likely to bear a star of more than a few tens of solar masses 
which would have accumulated on a timescale of ~ 10 5 years. 
However, the velocity structure around the protostars, and 
thus the instantaneous accretion rate, varies greatly from 
object to object, and this suggests that the final masses of 
the first stars may not lie in a narrow range. However, it 
is unclear whether the differences in the instantaneous ac- 
cretion rates will be directly proportional to the final mass 
of the stars because the velocity structure of the gas in the 
protostellar environment changes very rapidly with time. 

The fate of the star-forming clouds resolved in our sim- 
ulations is unclear. In order to reach stellar densities, the 
gas clouds must still collapse by more than four orders of 
magnitude in radius and this requires efficient transport of 
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angular momentum. This is a complicated problem and the 
various morphologies of the star-forming clouds we find sug- 
gest that its solution may involve more than a single physical 
process. One way to transport angular momentum outwards 
and facilitate further efficient accretion is through gravita- 
tional torques generated by spiral density fluctuations. Only 
3 of our 8 objects exhibit spiral structure in the central re- 
gions at the time when we halt our simulations and it is 
unclear if the other objects will also develop spiral structure 
later on. A further and critical complicating factor is the 
uncertain effect of feedback generated by the protostar on 
the dynamics of the accreting gas. 

Because of all of the complications just discussed, the 
masses of the first stars remain uncertain. A greater level 
of complexity than is possible with the current generation 
of simulations will be required to answer this fundamental 
question. These simulations will need to be able to model, in 
a self- consistent way, the dynamics of the accreting gas, the 
evolution of the protostar and the various radiative, chemical 
and mechanical feedback effects that are likely to be at play 
in the formation of the first generation of stars. 
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APPENDIX A: EVOLUTION OF VELOCITY 
PROFILES OF THE R5 OBJECT 

Figure Al is the evolution of velocity profiles of the R5 ob- 
ject. 
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Figure Al. Radially averaged profiles of velocity and angular 
momentum related physical quantities as a function of enclosed 
gas mass for 5 output times of the R5 object. The lines correspond 
to the same times as in Figure 3. Panel (a): Negative tangential 
velocity profiles. The dashed lines in this panel and panel (b) refer 
to negative sound speed. Panel (b): Radial velocity profiles. Panel 
(c): The ratio of free-fall time to cooling time. Panel (d): The ratio 
of tangential velocity to the required Keplerian velocity v Kepler — 
(GM/r) 1 / 2 . Panel (e): The orientation of angular momentum. 
Note that only the angle cf) is shown. The vertical doted lines in 
all panels indicate the enclosed gas mass within the virial radius. 
Note that only the last four outputs in panel (d) and (e) are 
shown, for clarity. 
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Figure Bl. The same plots as in Figure 7 but for a different set of 
5 simulations taken from our suite of runs: R5wt, Zlw, Zlwt, Z2 
and Z2w. Panel (a): Negative tangential velocity profiles. Panel 
(b): Radial velocity profiles. Panel (c): The ratio of the free-fall 
time to the cooling time. Panel (d): The ratio of the tangential 
velocity to the Keplerian velocity v Kepler — (GM/r) 1 / 2 . Panel 
(e): The orientation of the angular momentum profiles. Note that 
only the angle cf) is shown. The lines are the same as in Figure 6. 



APPENDIX B: VELOCITY PROFILES OF 5 
OBJECTS 

Figure Bl is the same plot as Figure 7 but are for other 5 
simulations. 
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